home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATH / FSULTRA1.ZIP / ULTRACOD.INC < prev    next >
Text File  |  1992-06-18  |  11KB  |  451 lines
  1. comment ! 
  2. FSU - ULTRA    The greatest random number generator that ever was
  3.         or ever will be.  Way beyond Super-Duper.
  4.         (Just kidding, but we think its a good one.)
  5.  
  6. Authors:    Arif Zaman (arif@stat.fsu.edu) and
  7.         George Marsaglia (geo@stat.fsu.edu).
  8.  
  9. Date:        27 May 1992
  10.  
  11. Version:    1.05
  12.  
  13. Copyright:    To obtain permission to incorporate this program into
  14.         any commercial product, please contact the authors at
  15.         the e-mail address given above or at
  16.  
  17.         Department of Statistics and
  18.         Supercomputer Computations Research Institute
  19.         Florida State University
  20.         Tallahassee, FL 32306.
  21.  
  22. See Also:    README        for a brief description
  23.         ULTRA.DOC    for a detailed description
  24.  
  25. -----------------------------------------------------------------------
  26. ! 
  27. ;
  28. ; This file defines the CODE segment to be included in a main file
  29. ; see ultra.doc for details.
  30. ;
  31. cong  macro     ; Congruential generator ==============================
  32. ;
  33. ; argument in dx:ax
  34. ; result   in dx:ax = low32bits of 69069 * arg
  35. ; clobbered   bx
  36. ;
  37. ; We need the lower 32 bits of 69069 * dx:ax
  38. ; Since we only have 16 bit multiplys with 32 bit answers,
  39. ; we use long multiplication. Let (a,b) represent a*2^16 + b.
  40. ;
  41. ;            69069 =                      (  1, 3533)
  42. ;            dx:ax =                 x    ( hi, lo  )
  43. ;                                 -------------------
  44. ;                     3533 * lo =         ( a1, a2  )
  45. ;                     3533 * hi =     ( b1, b2 )
  46. ;                        1 * lo =     (  0, lo )
  47. ;                                 -------------------
  48. ;            low 32 bits of answer = (a1+b2+lo, a2  )         
  49. ;
  50.       mov  bx,ax
  51.       mov  ax,3533
  52.       mul  dx                     ; dx:ax = b1,b2 = 3533*hi
  53.       xchg bx,ax
  54.       add  bx,ax                  ; bx = b2+lo;     ax = lo;
  55.       mov  dx,3533
  56.       mul  dx                     ; dx:ax = a1,a2 = 3533*lo
  57.       add  dx,bx                  ; dx = a1+b2+lo
  58.       endm
  59.  
  60. shreg macro     ; Shift register sequence =============================
  61. ;
  62. ;  argument is dx:ax
  63. ;  result   in dx:ax
  64. ;  clobbered   bx
  65. ;
  66. ;  consider the ax register as one sign bit `as' and the rest `a..r'
  67. ;  similarly decompose dx. The operation we want to do is:
  68. ;
  69. ;          x = x xor (x shr 15);
  70. ;          x = x xor (x shl 17);
  71. ;
  72. ;  Which can be shown in terms of the decomposed registers as:
  73. ;
  74. ;     x        = ds, d..r, as, a..r
  75. ;     x shr 15 =       ds, d..r  as
  76. ;                ------------------
  77. ;     x        = ds, d..r, as, a..r
  78. ;     x shl 17 = a..r   0
  79. ;                ------------------
  80. ;
  81. ;  The above operation is equivalent to:
  82. ;
  83. ;           ax = ax xor [d..r,as]
  84. ;           dx = dx xor [a..r,ds]
  85. ;
  86.       mov  bh,ah
  87.       shl  bh,1                 ; mov as into carry
  88.       mov  bx,dx
  89.       rcl  bx,1                 ; bx = [d..r,as]
  90.       pushf                     ; save carry flag (=ds)
  91.       xor  ax,bx                ; ax = ax xor bx
  92.       popf                      ; restore the carry
  93.       mov  bx,ax
  94.       rcl  bx,1                 ; bx = [a..r,ds]
  95.       xor  dx,bx                ; dx = dx xor bx
  96.       endm
  97.  
  98. rinit proc        ; INITIALIZE SEED ARRAY =======================
  99.     Enter2arg
  100. ;
  101. ; On exit:
  102. ;    The SWBseed array is set using sign bits from the McGill
  103. ;    Super-Duper generator, which is a mix of a congruential and
  104. ;    a shift-register sequence. Flags and counters are reset.
  105. ;
  106. ; Registers Clobbered:  AX,BX,CX,DX,SI,DI.
  107. ;
  108. ; Algorithm:
  109. ;    Starting from lsm of x[0] to the msb of x[N-1], each bit is
  110. ;    formed using the sign bit of the xor of a congruential generator
  111. ;    with seed ConX and a shift register sequence with seed ShrX.
  112. ;
  113. ; Register usage
  114. ;   cx contains count for two loops
  115. ;      ch counts the outer loop (for each byte of the array x)
  116. ;      cl counts for each bit of x[i]
  117. ;   bl contains the bits until eight of them are gathered
  118. ;   es:di point to where x[i] 
  119. ;   dx:ax are used for computations
  120. ;
  121.     mov ch,N*4
  122.     mov di,offset swbseed  ; do loop for x[0], ... , x[4*n] (bytes)
  123. nextbyte:
  124.     mov bl,0        ; 8 bits of bl are made bit by bit
  125.     mov cl,8        
  126. nextbit:    push bx                     ; compute the cong and the sh-reg
  127.         mov  ax,conglo              ; generators in place
  128.         mov  dx,conghi
  129.         cong
  130.         mov  conglo,ax
  131.         mov  conghi,dx
  132.         mov  ax,shrglo
  133.         mov  dx,shrghi
  134.         shreg
  135.         mov  shrglo,ax
  136.         mov  shrghi,dx
  137.         pop  bx 
  138.  
  139.         xor  dh,byte ptr conghi+1   ; xor the top byte of arg1 and arg2
  140.         rcl  dh,1           ; move the sign bit into the carry
  141.         rcr  bl,1           ; shift it into bl
  142.     dec cl
  143.     jnz nextbit
  144.     mov al,bl
  145.     stosb               ; Store the answer in swbb[i]
  146.     dec ch
  147.     jnz nextbyte
  148. ;
  149. ; reset all counters, flags etc. and return.
  150. ;
  151.     xor bx,bx
  152.     mov swb32.c,bx
  153.     mov swb16.c,bx
  154.     mov swb8.c,bx
  155.     mov swb1.c,bx
  156.     mov flags,0
  157.     Exit2arg
  158. rinit    endp
  159.  
  160. SWBfill proc near    ; SUBTRACT-WITH-BORROW ========================
  161. ;
  162. ;  On Entry:
  163. ;     DS should point to the data segment.
  164. ;     BX contains offset to swb??.x the data array for the answer
  165. ;
  166. ;  On Exit:
  167. ;     SWBseed contains all new values computed by the SWB generator.
  168. ;     carry   contains the carry flag from the last subtraction.
  169. ;     The array pointed to by [BX] contains the seed xor'ed with
  170. ;    a congruential sequence.
  171. ;
  172. ;  Registers Clobbered: AX,BX,CX,DX,SI,DI
  173. ;
  174. ;  Algorithm:
  175. ;     The following subtractions are performed from right to left.
  176. ;     The carry propagates as though these were two long numbers,
  177. ;     with each x[i] being a `digit' in base 2^32.
  178. ;
  179. ;        x[12] ...  x[ 0]      x[36] ...  x[13]
  180. ;       -x[36] ... -x[24]     -x[23] ... -x[ 0]
  181. ;       ---------------------------------------
  182. ;        x[36] ...  x[24]      x[23] ...  x[ 0]
  183. ;
  184. ;     The x's also could be considered as pairs of base 16 digits,
  185. ;     so that x[i] is the pair y[2i+1]y[2i]. This allows us to use
  186. ;     only 16 bit subtractions with carry, perfectly suited for all
  187. ;     80x86 processors. The same idea could be extended for machines
  188. ;     with only eight bit, or even only 1 bit arithmetic.
  189. ;
  190.     EnterFill
  191.     push bx            ; Save the location for the results
  192.     mov ah,byte ptr flags   ; set carry flag to what it was the
  193.     sahf                    ; last time we exited this routine
  194.  
  195.     mov cx,24*2             ; will do first loop 48 times
  196.     mov si,offset(swbseed)+13*4; set up ds:si -> x[13]
  197.     mov di,offset(swbseed)     ; set up es:di -> x[0]
  198.  
  199. loop1:  ; On a 8086, the instructions `lodsw', `stosw' and `loop'
  200.     ; all change registers automatically as noted in parentheses
  201.  
  202.     lodsw                   ;    ax = y[i+26]     ( si = si+2 )
  203.     sbb ax,[di]             ;    ax = ax-y[i]-carry
  204.     stosw                   ;    y[i] = ax        ( di = di+2 )
  205.     loop loop1              ; loop for i=0..47    ( cx = cx-1 )
  206.  
  207.     mov cx,13*2             ; will do next loop 26 times
  208.     mov si,offset(swbseed)     ; set up ds:si -> y[0]
  209.                     ; es:di is already set up
  210. loop2:
  211.     lodsw                   ;    ax = y[i-48]     ( si = si+2 )
  212.     sbb ax,[di]             ;    ax = ax-y[i]-carry
  213.     stosw                   ;    x[i] = ax        ( di = di+2 )
  214.     loop loop2              ; loop for i=48..73   ( cx = cx-1 )
  215.  
  216.     lahf
  217.     mov byte ptr flags,ah   ; save carry flag for next time
  218. ;
  219. ; XOR the elements of x with a congruential generator and put the
  220. ; result in swbx
  221. ;
  222.     mov  si,offset(swbseed)    ; the source is the seed
  223.     pop  di            ; the destination is one of swb??.x
  224.     mov  [di-4],di        ; set swb??.p = & swb??.x[0]
  225. IFNDEF fast
  226.     mov  cx,N
  227.     mov  ax,conglo
  228.     mov  dx,conghi
  229. loopc:  cong
  230.     mov  bx,ax
  231.     lodsw
  232.     xor  ax,bx
  233.     stosw
  234.     lodsw
  235.     xor  ax,dx
  236.     stosw
  237.     mov  ax,bx
  238.     loop loopc
  239.     mov  conglo,ax
  240.     mov  conghi,dx
  241. ELSE
  242.     mov  cx,2*N
  243.     rep  movsw
  244. ENDIF
  245.     ExitFill
  246. SWBfill    endp
  247.  
  248. ; Random Number procedures ============================================
  249. ;
  250. CheckFill MACRO bits,bytes,count
  251.     local ok
  252.     dec    bits.c
  253.     jns    ok
  254.         mov   bx,offset(bits.x)
  255.         call  SWBfill
  256.         mov   bits.c,count-1
  257. ok:    mov    bx,bits.p
  258.     add    bits.p,bytes
  259. ENDM
  260.  
  261. i32bit proc
  262.     EnterProcedure
  263.     CheckFill swb32,4,N
  264.     mov    ax,[bx]
  265.     mov    dx,[bx+2]
  266.     DwordFn
  267.     ExitProcedure
  268. i32bit endp
  269.  
  270. i31bit proc
  271.     EnterProcedure
  272.     CheckFill swb32,4,N
  273.     mov    ax,[bx]
  274.     mov    dx,[bx+2]
  275.     and    dh,7Fh
  276.     DwordFn
  277.     ExitProcedure
  278. i31bit endp
  279.  
  280. i16bit proc
  281.     EnterProcedure
  282.     CheckFill swb16,2,2*N
  283.     mov    ax,[bx]
  284.     WordFn
  285.     ExitProcedure
  286. i16bit endp
  287.  
  288. i15bit proc
  289.     EnterProcedure
  290.     CheckFill swb16,2,2*N
  291.     mov    ax,[bx]
  292.     and    ah,7Fh
  293.     WordFn
  294.     ExitProcedure
  295. i15bit endp
  296.  
  297. i8bit proc
  298.     EnterProcedure
  299.     CheckFill swb8,1,4*N
  300.     mov    al,[bx]
  301.     ByteFn
  302.     ExitProcedure
  303. i8bit endp
  304.  
  305. i7bit proc
  306.     EnterProcedure
  307.     CheckFill swb8,1,4*N
  308.     mov    al,[bx]
  309.     and    al,7Fh
  310.     ByteFn
  311.     ExitProcedure
  312. i7bit endp
  313.  
  314. i1bit proc
  315.     EnterProcedure di
  316.     dec    swb1.c
  317.     jns    ok1
  318.         CheckFill  swb32,4,N    ; do an i32bit
  319.         mov   dx,[bx+2]
  320.         mov   bx,[bx]        ; with answer in dx:bx
  321.         mov   swb1.c,31
  322.         mov   di,offset(swb1.x)
  323.         mov   swb1.p,di
  324.         mov   ax,ds
  325.         mov   es,ax
  326.  
  327.         mov   cx,16
  328. stosb1:        mov   al,1
  329.         and   al,bl
  330.         stosb
  331.         shr   bx,1
  332.         loop  stosb1
  333.  
  334.         mov      cl,16
  335. stosb2:        mov   al,1
  336.         and   al,dl
  337.         stosb
  338.         shr   dx,1
  339.         loop  stosb2
  340. ok1:    mov    bx,swb1.p
  341.     inc    swb1.p
  342.     mov    al,[bx]
  343.     ByteFn
  344.     ExitProcedure di
  345. i1bit endp
  346.  
  347. uni proc
  348.     EnterProcedure
  349.     fild    neg31
  350.     CheckFill  swb32,4,N
  351.     and    byte ptr [bx+3],7Fh
  352.     jnz    oku
  353.         mov      ax,[bx]
  354.         mov   tmpw3,ax
  355.         mov      ax,[bx+2]
  356.         mov      tmpw4,ax
  357.         CheckFill  swb32,4,N
  358.         mov      ax,[bx]
  359.         mov   tmpw1,ax
  360.         mov      ax,[bx]
  361.         mov   tmpw2,ax
  362.         fild  tmpq
  363.         jmp   uxit
  364. oku:    fild    dword ptr [bx]
  365. uxit:    fscale            ; scale the integer by the exponent
  366.     fstp st(1)        ; dump the exponent from the fpp
  367.     RealFn
  368.     ExitProcedure
  369. uni endp
  370.  
  371. vni proc
  372.     EnterProcedure
  373.     fild    neg31
  374.     CheckFill  swb32,4,N
  375.     test    byte ptr [bx+3],0FFh
  376.     jnz    okv
  377.         mov      ax,[bx]
  378.         mov   tmpw3,ax
  379.         mov   ax,[bx+2]
  380.         mov   tmpw4,ax
  381.         CheckFill  swb32,4,N
  382.         mov   ax,[bx]
  383.         mov   tmpw1,ax
  384.         mov   ax,[bx+2]
  385.         mov   tmpw2,ax
  386.         fild  tmpq
  387.         jmp   vxit
  388. okv:    fild    dword ptr [bx]
  389. vxit:    fscale            ; scale the integer by the exponent
  390.     fstp    st(1)        ; dump the exponent from the fpp
  391.     RealFn
  392.     ExitProcedure
  393. vni endp
  394.  
  395. duni proc
  396.     EnterProcedure
  397.     fild    neg63
  398.     CheckFill  swb32,4,N
  399.     dec    swb32.c
  400.     jns    okdu
  401.         mov      ax,[bx]
  402.         mov   tmpw1,ax
  403.         mov      ax,[bx+2]
  404.         mov   tmpw2,ax
  405.         mov   bx,offset(swb32.x)
  406.         call  SWBfill
  407.         mov   swb32.c,N-1
  408.         mov   ax,word ptr swb32.x
  409.         mov   tmpw3,ax
  410.         mov   ax,word ptr swb32.x+2
  411.         mov   tmpw4,ax
  412.         and   byte ptr [tmpq+7],7Fh
  413.         fild  tmpq
  414.         jmp   duxit
  415. okdu:    and    byte ptr [bx+7],7Fh
  416.     fild    qword ptr [bx]
  417. duxit:    add    swb32.p,4
  418.     fscale
  419.     fstp    st(1)
  420.     DoubleFn
  421.     ExitProcedure
  422. duni endp
  423.  
  424. dvni proc
  425.     EnterProcedure
  426.     fild    neg63
  427.     CheckFill  swb32,4,N
  428.     dec    swb32.c
  429.     jns    okdv
  430.         mov   ax,[bx]
  431.         mov   tmpw1,ax
  432.         mov   ax,[bx+2]
  433.         mov   tmpw2,ax
  434.         mov   bx,offset(swb32.x)
  435.         call  SWBfill
  436.         mov   swb32.c,N-1
  437.         mov   ax,word ptr swb32.x
  438.         mov   tmpw3,ax
  439.         mov   ax,word ptr swb32.x+2
  440.         mov   tmpw4,ax
  441.         fild  tmpq
  442.         jmp   dvxit
  443. okdv:    fild    qword ptr [bx]
  444. dvxit:    add    swb32.p,4
  445.     fscale
  446.     fstp    st(1)
  447.     DoubleFn
  448.     ExitProcedure
  449. dvni endp
  450.  
  451.